function [beta,out] = TypeRegressionEstimatorSmoothed(ThetaE,X_E,ThetaL,X_L,Q,contract,fem,blk,hsp,...
                                              domE,threshold1,threshold2,threshold3,...
                                              beta0,wts,UseGlobal,EqualVariance,EqualCorrelation,Display,tol_solver,TestExclusion)

if isempty(wts)
    wts = ones(size(ThetaE)) ;
end
rng(round(pi*1000000)) ;
Nsims = 1000 ;
SimZ = [randn(Nsims,1) randn(Nsims,1)] ;


K_E     = length(X_E(1,:)) ; 
K_L     = length(X_L(1,:)) ; 

if nargin<21
    TestExclusion = [] ;
elseif sum(size(TestExclusion)== [(K_E + K_L) 1]) ~= 2
    error('TypeRegressionEstimator:BadInput','TestExclusion must be a column with length=length(ThetaE)+length(ThetaL).')
elseif sum(TestExclusion==1|TestExclusion==0)~=length(TestExclusion)
    error('TypeRegressionEstimator:BadInput','TestExclusion must contain only 1s and 0s.')
end


tol_constraints = 0.05 ;
% tol_solver      = 10^-8 ;

%%%%LINEAR CONSTRAINTS: SIGMAe and SIGMAl must both be strictly positive and RHOel must be strictly
%%%%between -1 and 1.
Aineq = [zeros(1,K_E+K_L),-1,0,0,0,zeros(1,8)] ; %%%%baseline stdevE for white/asian/other males is positive
bineq = -tol_constraints/5 ;
Aineq = [Aineq; zeros(1,K_E+K_L),-1,-1,0,0,zeros(1,8)] ; %%%%stdevE for white/asian/other females is positive
bineq = [bineq; -tol_constraints/5] ;
Aineq = [Aineq; zeros(1,K_E+K_L),-1,0,-1,0,zeros(1,8)] ; %%%%stdevE for black males is positive
bineq = [bineq; -tol_constraints/5] ;
Aineq = [Aineq; zeros(1,K_E+K_L),-1,-1,-1,0,zeros(1,8)] ; %%%%stdevE for black females is positive
bineq = [bineq; -tol_constraints/5] ;
Aineq = [Aineq; zeros(1,K_E+K_L),-1,0,0,-1,zeros(1,8)] ; %%%%stdevE for hispanic males is positive
bineq = [bineq; -tol_constraints/5] ;
Aineq = [Aineq; zeros(1,K_E+K_L),-1,-1,0,-1,zeros(1,8)] ; %%%%stdevE for hispanic females is positive
bineq = [bineq; -tol_constraints/5] ;
Aineq = [Aineq; zeros(1,K_E+K_L),zeros(1,4),-1,0,0,0,zeros(1,4)] ; %%%%baseline stdevL for white/asian/other males is positive
bineq = [bineq; -tol_constraints/5] ;
Aineq = [Aineq; zeros(1,K_E+K_L),zeros(1,4),-1,-1,0,0,zeros(1,4)] ; %%%%stdevL for white/asian/other females is positive
bineq = [bineq; -tol_constraints/5] ;
Aineq = [Aineq; zeros(1,K_E+K_L),zeros(1,4),-1,0,-1,0,zeros(1,4)] ; %%%%stdevL for black males is positive
bineq = [bineq; -tol_constraints/5] ;
Aineq = [Aineq; zeros(1,K_E+K_L),zeros(1,4),-1,-1,-1,0,zeros(1,4)] ; %%%%stdevL for black females is positive
bineq = [bineq; -tol_constraints/5] ;
Aineq = [Aineq; zeros(1,K_E+K_L),zeros(1,4),-1,0,0,-1,zeros(1,4)] ; %%%%stdevL for hispanic males is positive
bineq = [bineq; -tol_constraints/5] ;
Aineq = [Aineq; zeros(1,K_E+K_L),zeros(1,4),-1,-1,0,-1,zeros(1,4)] ; %%%%stdevL for hispanic females is positive
bineq = [bineq; -tol_constraints/5] ;
Aineq = [Aineq; zeros(1,K_E+K_L+8), 1,0,0,0] ; %%%%baseline rho for white/asian/other males must be <1
bineq = [bineq; 1-tol_constraints] ;
Aineq = [Aineq; zeros(1,K_E+K_L+8), 1,1,0,0] ; %%%%rho for white/asian/other females must be <1
bineq = [bineq; 1-tol_constraints] ;
Aineq = [Aineq; zeros(1,K_E+K_L+8), 1,0,1,0] ; %%%%rho for black males must be <1
bineq = [bineq; 1-tol_constraints] ;
Aineq = [Aineq; zeros(1,K_E+K_L+8), 1,1,1,0] ; %%%%rho for black females must be <1
bineq = [bineq; 1-tol_constraints] ;
Aineq = [Aineq; zeros(1,K_E+K_L+8), 1,0,0,1] ; %%%%rho for hispanic males must be <1
bineq = [bineq; 1-tol_constraints] ;
Aineq = [Aineq; zeros(1,K_E+K_L+8), 1,1,0,1] ; %%%%rho for hispanic females must be <1
bineq = [bineq; 1-tol_constraints] ;
Aineq = [Aineq; zeros(1,K_E+K_L+8), -1,0,0,0] ; %%%%baseline rho for white/asian/other males must be >-1
bineq = [bineq; 1-tol_constraints] ;
Aineq = [Aineq; zeros(1,K_E+K_L+8), -1,-1,0,0] ; %%%%rho for white/asian/other females must be >-1
bineq = [bineq; 1-tol_constraints] ;
Aineq = [Aineq; zeros(1,K_E+K_L+8), -1,0,-1,0] ; %%%%rho for black males must be >-1
bineq = [bineq; 1-tol_constraints] ;
Aineq = [Aineq; zeros(1,K_E+K_L+8), -1,-1,-1,0] ; %%%%rho for black females must be >-1
bineq = [bineq; 1-tol_constraints] ;
Aineq = [Aineq; zeros(1,K_E+K_L+8), -1,0,0,-1] ; %%%%rho for hispanic males must be >-1
bineq = [bineq; 1-tol_constraints] ;
Aineq = [Aineq; zeros(1,K_E+K_L+8), -1,-1,0,-1] ; %%%%rho for hispanic females must be >-1 
bineq = [bineq; 1-tol_constraints] ;

lb = -inf*ones(1,K_E+K_L+12) ;
ub = inf*ones(1,K_E+K_L+12) ;

% EqualCorrelation = 1 ;
if EqualCorrelation==1
    Aeq = [zeros(1,K_E+K_L+8), 0,1,0,0] ;
    beq = 0 ;
    Aeq = [Aeq; zeros(1,K_E+K_L+8), 0,0,1,0] ;
    beq = [beq; 0] ;
    Aeq = [Aeq; zeros(1,K_E+K_L+8), 0,0,0,1] ;
    beq = [beq; 0] ;
    lb(1,K_E+K_L+8+1:K_E+K_L+8+4) = [tol_constraints-1,0,0,0] ;
    ub(1,K_E+K_L+8+1:K_E+K_L+8+4) = [1-tol_constraints,0,0,0] ;
else
    Aeq = [] ;
    beq = [] ;
end

% EqualVariance = 0 ;
if EqualVariance==1
    Aeq = [Aeq; zeros(1,K_E+K_L),0,1,0,0, zeros(1,8)] ;
    beq = [beq; 0] ;
    Aeq = [Aeq; zeros(1,K_E+K_L),0,0,1,0, zeros(1,8)] ;
    beq = [beq; 0] ;
    Aeq = [Aeq; zeros(1,K_E+K_L),0,0,0,1, zeros(1,8)] ;
    beq = [beq; 0] ;
    Aeq = [Aeq; zeros(1,K_E+K_L+4),0,1,0,0, zeros(1,4)] ;
    beq = [beq; 0] ;
    Aeq = [Aeq; zeros(1,K_E+K_L+4),0,0,1,0, zeros(1,4)] ;
    beq = [beq; 0] ;
    Aeq = [Aeq; zeros(1,K_E+K_L+4),0,0,0,1, zeros(1,4)] ;
    beq = [beq; 0] ;
    lb(1,K_E+K_L+1:K_E+K_L+4) = [tol_constraints,0,0,0] ;
    ub(1,K_E+K_L+2:K_E+K_L+4) = [0,0,0] ; 
    lb(1,K_E+K_L+5:K_E+K_L+8) = [tol_constraints,0,0,0] ;
    ub(1,K_E+K_L+6:K_E+K_L+8) = [0,0,0] ; 
end

%%%%THIS IS FOR TESTING JOINT EXCLUSION RESTRICTIONS VIA A LIKELIHOOD RATIO TEST:
if ~isempty(TestExclusion)
    beta0(TestExclusion==1) = 0 ;  
    Aeq = [Aeq; TestExclusion', zeros(1,12)] ;
    beq = [beq; 0] ;
end


fun = @objfun ;

disp('')
disp('Now running Tobit regression...')
% UseGlobal = 1 ;
if UseGlobal==0.1
    if Display==1
        options = optimoptions('fmincon','Algorithm','interior-point','FiniteDifferenceType','central',...
            'MaxIter',1000, 'MaxFunEvals',20000,...
            'TolCon',tol_solver, 'TolFun',tol_solver, 'TolX',tol_solver, 'Display','iter-detailed', 'UseParallel',true, 'HonorBounds',true) %#ok
    else
        options = optimoptions('fmincon','Algorithm','interior-point','FiniteDifferenceType','central',...
            'MaxIter',1000, 'MaxFunEvals',20000,...
            'TolCon',tol_solver, 'TolFun',tol_solver, 'TolX',tol_solver, 'Display','off', 'UseParallel',true, 'HonorBounds',true) ;
    end
    [beta,logLval,exitflag,output,lambda,grad,hess] = fmincon(fun, beta0, Aineq, bineq, Aeq, beq, lb, ub, [], options) ;
    out.lambda = lambda ;
    out.grad   = grad ; 
    out.hess   = hess ;
elseif UseGlobal==0.2
    if Display==1
        options = optimoptions('fmincon','Algorithm','sqp-legacy','FiniteDifferenceType','central', 'MaxIter',1000, 'MaxFunEvals',20000,...
            'TolCon',tol_solver, 'TolFun',tol_solver, 'TolX',tol_solver, 'Display','iter-detailed', 'UseParallel',true, 'HonorBounds',true) %#ok
    else
        options = optimoptions('fmincon','Algorithm','sqp-legacy','FiniteDifferenceType','central', 'MaxIter',1000, 'MaxFunEvals',20000,...
            'TolCon',tol_solver, 'TolFun',tol_solver, 'TolX',tol_solver, 'Display','off', 'UseParallel',true, 'HonorBounds',true) ;
    end
    [beta,logLval,exitflag,output,lambda,grad,hess] = fmincon(fun, beta0, Aineq, bineq, Aeq, beq, lb, ub, [], options) ;
    out.lambda = lambda ;
    out.grad   = grad ; 
    out.hess   = hess ;
elseif UseGlobal==0.3
    if Display==1
        options = optimoptions('fmincon','Algorithm','sqp','FiniteDifferenceType','central', 'MaxIter',1000, 'MaxFunEvals',20000,...
            'TolCon',tol_solver, 'TolFun',tol_solver, 'TolX',tol_solver, 'Display','iter-detailed', 'UseParallel',true, 'HonorBounds',true) %#ok
    else
        options = optimoptions('fmincon','Algorithm','sqp','FiniteDifferenceType','central', 'MaxIter',1000, 'MaxFunEvals',20000,...
            'TolCon',tol_solver, 'TolFun',tol_solver, 'TolX',tol_solver, 'Display','off', 'UseParallel',true, 'HonorBounds',true) ;
    end
    [beta,logLval,exitflag,output,lambda,grad,hess] = fmincon(fun, beta0, Aineq, bineq, Aeq, beq, lb, ub, [], options) ;
    out.lambda = lambda ;
    out.grad   = grad ; 
    out.hess   = hess ;
elseif UseGlobal==0.4
    if Display==1
        options = optimoptions('fmincon','Algorithm','active-set','FiniteDifferenceType','central', 'MaxIter',1000, 'MaxFunEvals',20000,...
            'TolCon',tol_solver, 'TolFun',tol_solver, 'TolX',tol_solver, 'Display','iter-detailed', 'UseParallel',true, 'HonorBounds',true) %#ok
    else
        options = optimoptions('fmincon','Algorithm','active-set','FiniteDifferenceType','central', 'MaxIter',1000, 'MaxFunEvals',20000,...
            'TolCon',tol_solver, 'TolFun',tol_solver, 'TolX',tol_solver, 'Display','off', 'UseParallel',true, 'HonorBounds',true) ;
    end
    [beta,logLval,exitflag,output,lambda,grad,hess] = fmincon(fun, beta0, Aineq, bineq, Aeq, beq, lb, ub, [], options) ;
    out.lambda = lambda ;
    out.grad   = grad ; 
    out.hess   = hess ;
elseif UseGlobal==1
    if Display==1
        options = optimoptions('patternsearch','Algorithm','classic', 'MaxIter',5000, 'MaxFunctionEvaluations',800000,...
            'MeshTolerance',tol_solver,'InitialMeshSize',16, 'Display','iter', 'UseParallel',true, 'PlotFcn',@psplotbestf) %#ok
    else
        options = optimoptions('patternsearch','Algorithm','classic', 'MaxIter',5000, 'MaxFunctionEvaluations',800000,...
            'MeshTolerance',tol_solver,'InitialMeshSize',16, 'Display','final', 'UseParallel',true) ;
    end
    [beta,logLval,exitflag,output] = patternsearch(fun, beta0, Aineq, bineq, Aeq, beq, lb, ub, [], options) ;
elseif UseGlobal==2  
    if Display==1
        options = optimoptions('patternsearch','Algorithm','nups', 'MaxIter',5000, 'MaxFunctionEvaluations',800000,...
            'MeshTolerance',tol_solver,'InitialMeshSize',16, 'Display','iter', 'UseParallel',true, 'PlotFcn',@psplotbestf) %#ok
    else
        options = optimoptions('patternsearch','Algorithm','nups', 'MaxIter',5000, 'MaxFunctionEvaluations',800000,...
            'MeshTolerance',tol_solver,'InitialMeshSize',16, 'Display','final', 'UseParallel',true) ;
    end
    [beta,logLval,exitflag,output] = patternsearch(fun, beta0, Aineq, bineq, Aeq, beq, lb, ub, [], options) ;
elseif UseGlobal==3  
    if Display==1
        options = optimoptions('patternsearch','Algorithm','nups-gps', 'MaxIter',5000, 'MaxFunctionEvaluations',800000,...
            'MeshTolerance',tol_solver,'InitialMeshSize',16, 'Display','iter', 'UseParallel',true, 'PlotFcn',@psplotbestf) %#ok
    else
        options = optimoptions('patternsearch','Algorithm','nups-gps', 'MaxIter',5000, 'MaxFunctionEvaluations',800000,...
            'MeshTolerance',tol_solver,'InitialMeshSize',16, 'Display','final', 'UseParallel',true) ;
    end
    [beta,logLval,exitflag,output] = patternsearch(fun, beta0, Aineq, bineq, Aeq, beq, lb, ub, [], options) ;
elseif UseGlobal==4  
    if Display==1
        options = optimoptions('patternsearch','Algorithm','nups-mads', 'MaxIter',5000, 'MaxFunctionEvaluations',800000,...
            'MeshTolerance',tol_solver,'InitialMeshSize',16, 'Display','iter', 'UseParallel',true, 'PlotFcn',@psplotbestf) %#ok
    else
        options = optimoptions('patternsearch','Algorithm','nups-mads', 'MaxIter',5000, 'MaxFunctionEvaluations',800000,...
            'MeshTolerance',tol_solver,'InitialMeshSize',16, 'Display','final', 'UseParallel',true) ;
    end
    [beta,logLval,exitflag,output] = patternsearch(fun, beta0, Aineq, bineq, Aeq, beq, lb, ub, [], options) ;
end
      

if nargout==2
    out.logLval = logLval ;
    out.exitflag = exitflag ;
    out.output   = output ;
end


    function logL = objfun(params)
        
        D0      = find(Q<2) ;  %%%%This collects all of the row indices for students with zero output
        D2      = find(Q>=2) ;  %%%find(Q>=2 & Q<80) ;  %%%%This collects all of the row indices for students with positive output
        lThetaE = log(ThetaE) ;  %%%%This is the log of ThetaE.  (ThetaE,ThetaL) is assumed to be a joint log-normal random variable.
        lThetaL = log(ThetaL) ;  %%%%This is the log of ThetaE.  (ThetaE,ThetaL) is assumed to be a joint log-normal random variable.
        %%%%The following compute the log of the selection thresholds:
        ldomE       = log(domE) ;
        lthreshold1 = log(threshold1) ;
        lthreshold2 = log(threshold2) ;
        lthreshold3 = log(threshold3) ;
        
        betaE = params(1:K_E) ;
        betaL = params(K_E+1:K_E+K_L) ;
        muE   = 0 ;
        muL   = 0 ;
        sigmaE     = params(K_E+K_L+1) ;
        sigmaE_f   = params(K_E+K_L+2) ;
        sigmaE_blk = params(K_E+K_L+3) ;
        sigmaE_hsp = params(K_E+K_L+4) ;
        SIGMAe     = sigmaE + sigmaE_f*fem + sigmaE_blk*blk + sigmaE_hsp*hsp ;
        sigmaL     = params(K_E+K_L+5) ;
        sigmaL_f   = params(K_E+K_L+6) ;
        sigmaL_blk = params(K_E+K_L+7) ;
        sigmaL_hsp = params(K_E+K_L+8) ;
        SIGMAl     = sigmaL + sigmaL_f*fem + sigmaL_blk*blk + sigmaL_hsp*hsp ;
        rhoEL      = params(K_E+K_L+9) ;
        rhoEL_f    = params(K_E+K_L+10) ;
        rhoEL_blk  = params(K_E+K_L+11) ;
        rhoEL_hsp  = params(K_E+K_L+12) ;
        RHOel      = rhoEL + rhoEL_f*fem + rhoEL_blk*blk + rhoEL_hsp*hsp ;
        
        logL = zeros(length(X_L(:,1)),1) ;  %%%%This will hold the contributions to the log-likelihood function, which we evaluate in three parts...
        
        %%%%PART 1: all rows where both ThetaE and ThetaL are known
        X2 = [(lThetaE(D2)-X_E(D2,:)*betaE) (lThetaL(D2)-X_L(D2,:)*betaL)] ;
        for ii=1:length(D2)
            S = [SIGMAe(D2(ii))^2, RHOel(D2(ii))*SIGMAe(D2(ii))*SIGMAl(D2(ii)); ...
                RHOel(D2(ii))*SIGMAe(D2(ii))*SIGMAl(D2(ii)), SIGMAl(D2(ii))^2] ;
            try
            logL(D2(ii)) = log( mvnpdf(X2(ii,:),[0,0],S) ) ;
            catch
                save('SNAFU.mat')
                error()
            end
        end
        
        
        %%%%PART 3: all rows where Q<2 (i.e., where (ThetaE,ThetaL) lives in a bounded region but
        %%%%neither is known.
        for ii=1:length(D0)
            M = [muE muL] ;
            S = [SIGMAe(D0(ii))^2, RHOel(D0(ii))*SIGMAe(D0(ii))*SIGMAl(D0(ii)); ...
                RHOel(D0(ii))*SIGMAe(D0(ii))*SIGMAl(D0(ii)), SIGMAl(D0(ii))^2] ; 
            
            %%%%The following code transforms a 2-D vector of iid standard normal values into
            %%%%general random variables with mean vector M and variance-covariance matrix S.
            %%%%STEP 1: compute the Cholesky factorization of the var-cov matrix:
            V = chol(S) ;  %%%%NOTE: S = V'*V
            %%%%STEP 2: here is how we convert all of the Z-values in SimZ to general Normal RVs
            SimNorm = ( V'*SimZ' + M'*ones(1,Nsims) )' ;
            %%%%STEP 3: add the regressors to the simulated residuals with the appropriate
            %%%%distribution from above and we have complete simulations from (ThetaE,ThetaL) space.
            SimNorm(:,1) = SimNorm(:,1) + X_E(D0(ii),:)*betaE ;
            SimNorm(:,2) = SimNorm(:,2) + X_L(D0(ii),:)*betaL ;
            % % % % % bandwidth = 1.06*2.978*std(AchievementData.ThetaL(~isnan(AchievementData.ThetaL)))*(sum(~isnan(AchievementData.ThetaL))^(-0.2)) ;
            bandwidth = 0.154173746772965 ;
            if contract(D0(ii))==1
                % % % % % outflag = 1*( SimNorm(:,2) > interp1(ldomE,lthreshold1,SimNorm(:,1),'linear','extrap') ) ;
                outflag = triweightCDF( (SimNorm(:,2) - interp1(ldomE,lthreshold1,SimNorm(:,1),'linear','extrap'))/bandwidth ) ;
                PrOut   = mean(outflag) ;
                if PrOut==1
                    PrOut = 1 - 1/(Nsims+1) ;
                elseif PrOut==0
                    PrOut = 1/(Nsims+1) ;
                end
                logL(D0(ii)) = log(PrOut) ;
            elseif contract(D0(ii))==2
                % % % % % outflag = 1*( SimNorm(:,2) > interp1(ldomE,lthreshold2,SimNorm(:,1),'linear','extrap') ) ;
                outflag = triweightCDF( (SimNorm(:,2) - interp1(ldomE,lthreshold2,SimNorm(:,1),'linear','extrap'))/bandwidth ) ;
                PrOut   = mean(outflag) ;
                if PrOut==1
                    PrOut = 1 - 1/(Nsims+1) ;
                elseif PrOut==0
                    PrOut = 1/(Nsims+1) ;
                end
                logL(D0(ii)) = log(PrOut) ;
            elseif contract(D0(ii))==3
                % % % % % outflag = 1*( SimNorm(:,2) > interp1(ldomE,lthreshold3,SimNorm(:,1),'linear','extrap') ) ;
                outflag = triweightCDF( (SimNorm(:,2) - interp1(ldomE,lthreshold3,SimNorm(:,1),'linear','extrap'))/bandwidth ) ;
                PrOut   = mean(outflag) ;
                if PrOut==1
                    PrOut = 1 - 1/(Nsims+1) ;
                elseif PrOut==0
                    PrOut = 1/(Nsims+1) ;
                end
                logL(D0(ii)) = log(PrOut) ;
            end
        end
        
        logL = -sum(wts.*logL) ; 
    end

    function F = triweightCDF(x) 
        F           = zeros(size(x)) ;
        F(x>1)      = 1 ;
        betwindx    = (x>-1 & x<=1) ;
        F(betwindx) = (35/32)*(x(betwindx) - x(betwindx).^3 + (3/5)*x(betwindx).^5 - (1/7)*x(betwindx).^7) + (1/2) ; 
    end


end



